/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		fem.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-04 15:40
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:27
*
*   Description:
*
=============================================================================*/
#ifndef FEM_H
#define FEM_H

#include "geometry/mesh/build_mesh.hpp"
#include "model/solve/assemble_equation.hpp"

namespace geoxfem {

namespace model {
struct fem {};

template <typename CouplePolicy>
struct assemble_equation<fem, CouplePolicy> { /*{{{*/
  template <typename Mesh, typename DOFHandler, typename TimeStep>
  static void apply(Mesh& mesh, DOFHandler& dof_handler, TimeStep& tstep) {
    size_t dimension = traits::dimension<Mesh>::value;
    auto&  kay_mat   = tstep.getSparseKay();
    auto&  rhs       = tstep.getRhs();
    auto&  answer    = tstep.getAnswer();
    GEOXFEM_LOG_SEV(info) << "assembling global kay matrix and rhs...";
    /* 1. reset and malloc arrays and matrix according to the number of Dofs */
    auto neq = static_cast<Eigen::Index>(dof_handler.getNEQ());
    kay_mat.resize(neq, neq);
    kay_mat.data().squeeze();
    rhs.setZero(neq);
    answer.setZero(neq);

    /* 2. add external load from Dofs into RHS vector(neumann boundary) */
    for(auto& dof : dofsIn(dof_handler))
      if(dof->ifNeumann())
        rhs[dof->getIndex()] += dof->getNeumannValue();
        /* 3.1 assemble global kay matrix */
        /* 3.2 deal with dirichlet   boundary(modify kay matrix and RHS) */
        /* 3.3 form iteration step load into rhs from internal state */
#pragma omp parallel
    {
      Eigen::SparseMatrix<_real> kay_mat_private;
      kay_mat_private.resize(neq, neq);
      VectorXR rhs_private;
      rhs_private.setZero(neq);
#pragma omp for nowait
      for(auto it = mesh.getShapeBegin(dimension);
          it != mesh.getShapeEnd(dimension); ++it)
        assemble_element<CouplePolicy>::apply(**it, tstep, kay_mat_private,
                                              rhs_private);
#pragma omp critical
      {
        // Here we sum the results of each thread one at a time
        kay_mat += kay_mat_private;
        rhs += rhs_private;
      }
    }
    /* GEOXFEM_ASSERT(!kay_mat.hasNaN()); */
    GEOXFEM_ASSERT(!rhs.hasNaN());
    return;
  }
}; /*}}}*/

}  // namespace model

namespace geometry {

template <>
struct build_mesh_impl<model::fem> { /*{{{*/
  using ptree = boost::property_tree::ptree;

  template <typename Model, typename Mesh>
  static void apply(const ptree& root, Model& model, Mesh& mesh) {
    using DimensionType = typename Model::DimensionType;
    using CouplePolicy  = typename Model::CouplePolicy;
    using Shape         = typename Model::Shape;
    using point_type    = typename traits::point_type<DimensionType>::type;
    using Material      = material::MaterialInterface<DimensionType>;
    {
      GEOXFEM_LOG_SEV(info) << "Building nodes from source file...";
      /* 1. build nodes, iterator all nodes record from xml source file */
      for(auto& node_data : root.get_child("nodes")) {
        point_type point = model.append(point_type(node_data.second));
        std::unique_ptr<Shape> spVertex =
          Shape::getFactory().create("Vertex", {point});
        set_nodal_dofs<CouplePolicy>::apply(spVertex->toVertex());
        GEOXFEM_LOG_SEV(trace) << "Building " << *spVertex;
        model.append(spVertex, &mesh);
      }
    }
    {
      GEOXFEM_LOG_SEV(info) << "Building elements from source file...";
      /* 2. build element, iterator all elements record from xml source file */
      for(auto& elem_data : root.get_child("elements")) {
        /* 2.1 assemble shape from input data */
        std::string type = elem_data.second.template get<std::string>("type");
        std::vector<point_type> node_list;
        for(auto& v : elem_data.second.get_child("nodes")) {
          size_t     point_id = std::stoul(v.second.data());
          point_type point    = model.getPoint(point_id).shallowClone();
          node_list.emplace_back(point);
        }
        std::unique_ptr<Shape> spShape =
          Shape::getFactory().create(type, node_list);
        size_t id = elem_data.second.template get<size_t>("id");
        spShape->setIndex(id);
        GEOXFEM_LOG_SEV(trace) << "Building " << *spShape;
        int mat_id = elem_data.second.template get<int>("mat", -1);
        if(mat_id > -1) {
          const Material* mat = model.getMaterial(mat_id);
          spShape->setMaterial(mat);
        }
        /* 2.2 add new shape into model and given mesh */
        auto shp = model.append(spShape, &mesh);
        /* 2.3 accumulate nodal weight */
        shp->setNodalWeight();
      }
    }
    return;
  }
}; /*}}}*/

}  // namespace geometry
}  // namespace geoxfem

#endif /* FEM_H */
